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We investigate finite size effects in quantum quenches on the basis of simple energetic arguments. Distin- 
guishing between the low-energy part of the excitation spectrum, below a microscopic energy-scale, and the 
high-energy regime enables one to define a crossover number of particles that is shown to diverge in the small 
quench limit. Another crossover number is proposed based on the fidelity between the initial and final ground- 
states. Both criteria can be computed using ground-state techniques that work for larger system sizes than full 
spectrum diagonalization. As examples, two models are studied: one with free bosons in an harmonic trap which 
frequency is quenched, and the one-dimensional Bose-Hubbard model, that is known to be non-integrable and 
for which recent studies have uncovered remarkable non-equilibrium behaviors. The diagonal weights of the 
time-averaged density-matrix are computed and observables obtained from this diagonal ensemble are compared 
with the ones from statistical ensembles. It is argued that the "thermalized" regime of the Bose-Hubbard model, 
previously observed in the small quench regime, experiences strong finite size effects that render difficult a thor- 
ough comparison with statistical ensembles. In addition, we show that the non-thermalized regime, emerging 
on finite size systems and for large interaction quenches, is not related to the existence of an equilibrium quan- 
tum critical point but to the high energy structure of the energy spectrum in the atomic limit. Its features are 
reminiscent of the quench from the non-interacting limit to the atomic limit. 
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The study of the non-equilibrium evolution of closed quan- 
tum many-body systems has been triggered by the recent pro- 
gresses in cold atoms experiments in which atoms are hardly 
coupled to the environment (Ql, 0] ■ Furthermore, microscopic 
parameters of the Hamiltonian governing the dynamics can 
be controlled at will and changed on microscopic timescales. 
In this context, the question of the unitary evolution of an iso- 
lated quantum system after a sudden change of one parameter, 
the so-called quantum quench, has attracted a lot of interest in 
both the experimental and theoretical communities [jst]. Many 
different questions are raised by such a set-up, among which 
are the relaxation of observables |4- 17], the question of ther- 
malization j 1 li 1 1 8l434ll . the existence of a subsystem steady- 
state lIssUssll . and the propagation of the entanglement |39- 
I43I1 . Beyond these academic concerns, practical applications 
of quenches have been proposed through the engineering of 
metastable states 1,44. ,45] and of an out-of-equilibrium super- 
solid state in a cold atoms set-up ll46ll . This paper is dedi- 
cated to the thermalization issue, but restricted to specific ex- 
amples and without claims on general results about the ther- 
malization mechanism. In this context, a quench can be un- 
derstood as a way to create an initial state that evolves through 
the dynamics of a given Hamiltonian. A common wisdom in 
classical mechanics is that the long-time evolution will for- 
get about the initial state and will explore all the accessible 
phase-space, provided the dynamics are chaotic. Then, er- 
godicity allows for the use of statistical ensembles in place 
of time-averaging. For a closed quantum system, as the evo- 
lution is unitary and the spectrum discrete, long-time recur- 
rences occur and the contribution of the eigenstates involved 
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in the dynamics is fixed by the initial state. For large enough 
systems, a quantum ergodic theorem was proposed |47], sup- 
porting the emergence of the microcanonical ensemble which 
is the usual statistical ensemble for an isolated system. This 
approach aims at showing that time-averaged density-matrix p 
(see below for the definition) is macroscopically equivalent to 
the microcanonical ensemble. In a quantum quench, the initial 
state is not a "typical" state of a given energy but usually is the 
ground-state of the same Hamiltonian with different parame- 
ters. Consequently, the quench amplitude, or how much do 
we change the Hamiltonian, is here another relevant quantity. 
Another way to regard a quench can be as a perturbation of 
the initial state and one may wonder whether the long-time re- 
sponse is sensitive to the initial state. Furthermore, numerical 
tools and experiments on closed systems cannot easily reach 
a large number of particles so finite size effects can be im- 
portant in the interpretation of the observed phenomena. This 
paper suggests possible approaches to the question of these 
finite size effects after a quantum quench, and a possible in- 
terpretation of the observations made on a particular model: 
the one-dimensional Bose-Hubbard model. The other model, 
consisting of free bosons in an harmonic trap, offers another 
example of finite size effects and remarkable behaviors. Some 
of the features of the two models are surprisingly connected. 

The central object governing the long-time physics after a 
quantum quench is the time-averaged density-matrix p that 
predicts the time-averaged expectation values of any observ- 
able. This density-matrix has also connections to the heat or 
work done on a system 1.4.8.-5 1 ] . The weights of this "diag- 
onal ensemble" are difficult to compute for large systems as 
one needs to fully diagonalize the Hamiltonian so one unfortu- 
nately has to work with small systems (Hilbert spaces). Other 
methods have been used to tackle the physics of quenches. 
For instance, "Ab-initio" numerics have been used on both in- 
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tegrable and non-integrable models ifTH [Tl [T9ll27i [ssIIstI 
l38l I52I - I55I1 . Numerical methods like time-dependent density- 
matrix renormalization group (tDMRG) |56-58] can be used 
to compute the time-evolution of the wave-function but the in- 
terpretation is restricted to observables and to a finite window 
of time, and cannot give access the these weights. Exact re- 
sults on integrable models fl HI [H [H H El HiMH have 
the advantage to treat in a non-perturbative way large systems, 
but on the other hand, it is not surprising that they do not 
always thermalize due to the extensive number of conserved 
quantities. Luttinger liquid theory, which describes the low- 
energy physics of one-dimensional models in terms of free 
bosonic fields (thus an integrable theory), has been used to 
compute the time-evolution of the observables ll7l-l9l l62[l63ll . 
Quantum chaos methods have also helped studying the time- 
evolution of the Bose-Hubbard model [64-66]. Some stud- 
ies focused on the relation between fidelity and on the energy 
distribution lf54ll67ll . All these methods suffer from approxi- 
mations and/or finite size effects and it is sometimes hard to 
determine what is an artifact or not. 

Some of the results from numerical simulations seems to 
be contradictory IH [Ull^ IH, [H HI but were carried 
out on different models with different range of parameters, 
and not necessarily starting from the ground-state |27] of a 
simply related Hamiltonian. Performing a quantum quench 
amounts to projecting an initial state onto the energy spectrum 
of the final Hamiltonian, corresponding to a certain distribu- 
tion of energy p{E). In the thermodynamical limit, a global 
quench is expected to drive the mean energy to the bulk of 
the energy spectrum since the perturbing operator is exten- 
sive. In this high energy domain, semi-classical physics and 
random-matrix theory arguments are expected to work and 
make expectation values hardly depend on the energy^ithin 
a window given by the energy fluctuations) |18, 23, I21II : 
thermalization can occur in the sense that the energy dis- 
tribution obtained from the quench gives the same averages 
for the observables as the microcanonical ensemble. This 
so-called "eigenstates thermaUzation hypothesis" (ETH) has 
been tested numerically ifTol IztI [32I [33:1 for given models 
(typically fermionic and hard-core bosonic models) and some 
given set of parameters. No memory of the initial state 
(for a given mean energy) is thus found on simple observ- 
ables. These results agree well with the previous findings of 
Ref. 121 on a similar model. Having in mind this qualitative 
argument, the results of Ref. 11 on the non-integrable one- 
dimensional (ID) Bose-Hubbard model (BHM) look rather 
counter-intuitive: for small quenches, a thermalized regime 
was found in the sense that two independent observables com- 
puted within a (grand)-canonical ensemble (and not micro- 
canonical) and from time-evolution gave the same results. On 
the contrary, a mean-field treatment of the ID BHM inter- 
preted in the framework of chaos theory (6^] supports non- 
thermalization below an interaction threshold and thermaliza- 
tion above (mean-field theory is however known to fail for this 
strongly-correlated model so the results are not under control). 
The findings of Ref. 1 1 were later supported by the calculation 
of the diagonal ensemble distributions which looked like an 
approximate Boltzmann law [541 in the small quench regime. 



Surprisingly, for large quenches, a non-thermalized regime 
was found in Ref. 11 in which the correlations bear a strong 
memory of the initial state (in the sense that they are closer to 
the ones in the initial state than to the thermalized ones). This 
non-equilibrium behavior was attributed to the very peculiar 
shape of the diagonal ensemble in this regime 15411 . An im- 
portant step towards the understanding of the non-thermalized 
regime on finite size systems was made very recently [551 by 
giving numerical evidences on the ID BHM that ETH does 
not apply for large quenches in finite systems and suggesting 
a general framework in terms of rare events contributing to the 
distribution, providing a refined version of the ETH. 

As integrability is often one of the ingredients that play 
a role in the physics of quenches, we briefly recall that, for 
ID quantum many -body models, integrability can be well de- 
fined for a class of models which have the property of scatter- 
ing without diffraction ll68ll . This has two consequences that 
are in relation with the question of thermalization: the mo- 
menta of the particles do not redistribute [681 (a process which 
is believed to be essential to get the thermalized momentum 
distribution), and there is an extensive number of conserved 
quantities that separate the eigenstates in many sectors, con- 
straining the time evolution. In the context of nuclear physics, 
random-matrix theory has been proposed to describe the sta- 
tistical features of the bulk of the spectrum and it is commonly 
conjectured that non-integrable quantum many-body or clas- 
sically chaotic models display universal level statistics ll69ll . 
Level statistics have been computed in a few many -body mod- 
els [701, supporting the conjecture, but these results are re- 
stricted to a few models and it cannot be excluded that diffrac- 
tive models could display non-universal level statistics. The 
Bose-Hubbard model is a bit peculiar in this sense: if one 
denotes by A^max the maximum number of bosons onsite, the 
model is non-diffractive only for A'^ax = 1 |711. In addi- 
tion, if U is the interaction strength, U = is an integrable 
point (the atomic limit J = is as well exactly solvable). 
Level statistics and delocalization properties of the eigenstates 
have shown iFtTI I72I1 that the BHM display features of quan- 
tum chaotic systems for non-zero U (and larger A^^ax)- 

The first goal of this paper is to discuss the crossover from 
small to large quench amplitude regimes on the basis of en- 
ergetic and static fidelity arguments, and to evaluate the finite 
size effects that are associated to this crossover We then turn 
to a detailed discussion of the diagonal ensemble and the ver- 
ification of the ETH in the BHM, complementary to what has 
been done in Refs. |54J and |55). We show that the observed 
Boltzmann-like regime is spoiled by strong finite size effects 
that prevent both an accurate definition of an effective temper- 
ature and the comparison with the microcanonical ensemble. 
In the large quench limit, we explain in details that the break- 
down of the ETH is actually related to the "integrable" quench 
limit [/i = —>[// = 00. Thus, non-thermalization in the 
ID BHM is, on finite systems, reminiscent of the atomic limit. 
While the U = limit of the Bose-Hubbard model is trivially 
integrable as a free boson model, the infinite U (or atomic) 
limit is a bit particular: for very large U and focusing on the 
low-energy part of the spectrum, the model is effectively iden- 
tical to an integrable ID hard-core bosons model (A^ax = !)■ 
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However, we will see that, to understand the large-i7 limit 
of the quench, we will have to consider the whole excitation 
spectrum and not only the low-energy part. This result can be 
qualitatively and partially connected to the effect of the prox- 
imity to integrable points in quantum quenches, studied very 
recently in fermionic and hard-core bosonic models ll32[[33ll . 
in the sense that the observed non-thermalized regime on finite 
systems is connected to a particular limit in which the model 
has high degeneracies. Throughout the paper, we also give 
a simple but interesting example of a quench in a toy model 
consisting of free bosons confined in an harmonic trap. The 
motivation for it is that it surprisingly shares some qualitative 
features with the ID BHM and that it allows for analytical 
calculations on some properties of the diagonal ensemble dis- 
tribution. This model also corresponds to a standard exper- 
imental setup (so for the BHM) although interactions would 
have to be taken into account for a realistic comparison. 

The paper is organized as follows: we first review in Sec. U 
the definitions of the time-averaged density-matrix, the ETH 
and the computation of the diagonal weights for the two mod- 
els under study. In Sec.|ll] we suggest two kinds of crossover 
number of particles to distinguish the small and large quench 
regimes. Lastly, we discuss in Sec.|III]the fate of the ETH in 
the ID BHM and on small finite size systems. 



MODELS AND COMPUTATION OF THE WEIGHTS OF 
THE DIAGONAL ENSEMBLE 

A. The time-averaged density-matrix and the "eigenstate 
thermaUzation hypothesis" 



As discussed in recent papers 11271 132 |, I33L 154 155L |60|] , the 

time-averaged expectation values of any observable are gov- 
erned by the time-averaged density-matrix p, which is diag- 
onal in the final Hamiltonian eigenstate basis, provided the 
spectrum is non-degenerate. From now on, we only consider 
finite size systems that have a discrete spectrum. This leads 
to the so-called "diagonal ensemble" that has weights fully 
determined by the overlaps between the initial state IV'o.i) 
and eigenstates l^/'n./) of the final Hamiltonian H/. Usu- 
ally, \'4}Q.i) is the ground-state of the initial Hamiltonian "H; 
and we assume in the following that we start from this zero- 
temperature pure state. We also consider that the final Hamil- 
tonian takes the form 



(1) 



where A (that has the dimension of an energy) is called 
the quench amplitude, and T-Li is the dimensionless perturb- 
ing operator. Working on a global quantum quench means 
that "Hi is assumed to be an extensive operator that scales 
with the number of particles N . The time-averaged density- 
matrix is defined hy p ~ limt^oo \ Jq IV'('S)) (V'(s)l ds with 
\ip{t)) — e^*^^* IV'o.j)- It is important to realize that the in- 
finite time limit is taken before the thermodynamical limit. 
If the spectrum has exact degeneracies, the time-averaged 



density-matrix reads: 

p = J2p^ (^"jI + i^-*'/) ^^-i'/i 

n d 

where n labels non-degenerate eigenstates of Hf and p„ = 
|('0n,/|'0o,i)P are the diagonal weights, d labels the ba- 
sis of the degenerate subspaces, and the vectors \ipd,f) = 
J2q^ f{ldj\i^o,i) \qdj) keep a memory of the initial phases 
of I'i/'o.i) with respect to the \qd.f)- In the situation where p is 
block-diagonal, in order to get time-averaged results for an ob- 
servable O which has off-diagonal matrix elements in the H f 
eigenstate basis, one would have to compute all the overlaps 
('?d,/|V'o,i) ™d {q'^ j\0\qd.f) and sum up the contributions of 
all a degenerate subspace. In the following, this would be the 
case only for the free boson model and we will actually only 
use observables that are diagonal because the dimensions of 
the degenerate sectors grows (roughly) exponentially with the 
number of bosons N. For the Bose-Hubbard model, one can 
check that the spectra are non-degenerate in each symmetry 
sector. 

For a generic non-integrable model, the "eigenstate ther- 
malization hypothesis" (ETH) has been surmised I20l - l23ll27ll . 
suggesting an explanation for thermalization in an isolated 
quantum system and a justification for the use of the micro- 
canonical ensemble. The ETH is supported by semi-classical 
and random-matrix theory arguments ifTsl uOi- 23ll. and was 
checked numerically on particular models 1119 , l27l [32i [33ll . 
The ETH boils down to the fact that, in a given small win- 
dow of energy, the diagonal observables 0„ = {lpnj\0\lpn,f) 
that contribute to the time-averaged expectation value O = 
Tr[pC'] = J2ja PnOn hardly depend on the eigenstate n (in 
short. On ~ O in a small energy window). Consequently, any 
distribution peaked around the mean energy, and one can show 
on general grounds that the relative width of the distribution 
scales to zero as 

TV- 1271] (although some slower scalings 
could occur ifs^l t. will give the same observables as the micro- 
canonical ensemble, therefore accounting for thermalization. 
For integrable models ll27l[32l[33ll . non-thermalization is ex- 
plained from the fact that observables fluctuate a lot within 
a given energy window, which may be associated with the 
extensive number of conserved quantities that exist in these 
models. A more subtle scenario for the breakdown of the ETH 
was recently proposed i55[l . in which some "rare" states have 
a significant contribution to the averaged observables. 



B. Free bosons in an harmonic trap 

We now describe how to get the diagonal weights for two 
particular models. Firstly, we consider a model of N non- 
interacting bosons initially confined in an harmonic trap of 
frequency and lying in the zero-temperature ground-state. 
The frequency is changed to at time t — Q. For this 
model, the quench amplitude is defined as A = ujf/uji — 1 
(taking uji as the unit of energy), according to the expression 
of the quench parameter in terms of the harmonic oscillator 
ladder operators. We start with the computation of the single- 
particle wave-function overlaps p„ since the results for the 
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many-body wave-function will be expressed as a function of 
them. The single-particle spectrum is non-degenerate and the 
single-particle eigenfunctions are: 



following Hamiltonian: 



1 _ x2 /X 

—====e '^B.A- 



with cr = ^/h/rnuj and H„ the Hermite polynomials. The 
single-particle excitation spectrum is split into the odd and 
even parity sectors and the overlaps are non-zero for even- 
parity wave-functions only. They read: 



P2n 



(2n)! VT+A / A 



22»(n!)2 l + A/2 VA + 2 



(3) 



for integer n. The many-body wave-function of an A^-bosons 
excited configuration {fij} — {uq,--- ,nm} of the final 
Hamiltonian T-l f (with highest occupied level to) is: 



no!ni! • • 



Nl 



with V the set of all permutations and rij the occupation of 
the single-particle orbital (j)j.f- Overlapping this state with the 
iV-bosons initial ground-state |0o.i, • • • , 0o.i) gives the many- 
body weights 



Nl 



n2\ 



(4) 



In this equation, all to's are even integers. The total energy of 
this excitation is E^n.y = hujf{2n2 + 4n4 + • • • + mn„i) + 

hjjfN/2 with the constraint Yl^lo = N. Eq. (|4|l is noth- 
ing but the multinomial distribution associated with the ele- 
mentary probabilities p„, and it is thus clear that it is normal- 
ized. We also see that formula Q is in general valid for a free 
boson model, starting from the condensed ground-state (and 
specifying the pm)- If one takes the single-particle Boltzmann 
factor for the p™, one recovers the many-body Boltzmann fac- 
tor for the configuration. Contrary to statistical ensemble dis- 
tributions, the weights do not show a simple dependence of 
the configuration energy. This quench is qualitatively similar 
to a Joule compression/expansion as the ID effective density 
n — Nll! suddenly changes. In fact, A = n//ni — lis related 
the ratio of the effective densities. Other examples of quantum 
mechanical treatments of the Joule expansion can be found in 
the literature ll73ll74ll . 

In order to get the distribution of the weights versus en- 
ergy, we resort to numerics: using a fixed number of low-lying 
even parity levels Ns, we scan all possible configurations of 
N bosons in these Ng levels iteratively up to roughly 62 x 10^ 
configurations (N = 18 and Ng = 22). The truncation er- 
ror associated with a finite Ng is checked by summing up the 
weights. 



C. The one-dimensional Bose-Hubbard model 

The Bose-Hubbard model in a one-dimensional lattice, 
known to be non-integrable for [/ 7^ 0, is described by the 



•H = - J 



with the operator creating a boson at site j and nj ~ h^-bj 
the local density. J is the kinetic energy scale while U is the 
magnitude of the onsite repulsion. In an optical lattice, the 
ratio U / J can be tuned by changing the depth of the lattice 
and using Feshbach resonance. When the density of bosons 
is fixed at n = 1 and U is increased, the zero-temperature 
equilibrium phase diagram of the model displays a quantum 
phase transition from a superfluid phase to a Mott insulating 
phase in which particles are localized on each site. The criti- 
cal point has been located at Uc — 3.3 J using numerics iFtsIi . 
The quenches are performed by changing the interaction pa- 
rameter Ui ^ Uf (we set J = 1 as the unit of energy in the 
following), so we have A = {Uf — Ui)/2, and the perturb- 
ing operator "Hi = J^j ~ 1) is diagonal. Numerically, 
one must fix a maximum onsite occupancy A'^ax and we take 
A'max = 4 unless stated otherwise. Exact diagonalization cal- 
culations are carried out using periodic boundary conditions 
and translational invariance. We denote byO<fc<L — 1 the 
total momentum symmetry sectors. The algorithm to get the 
ground-state and eigenstates of the Hamiltonian is a full diag- 
onalization scheme for sizes up to L = 10 at unitary filling. 
For some of the quantities, we use the Lanczos algorithm up 
to L = 15. In Ref. 54, the Lanczos algorithm has been pro- 
posed to compute the low-energy weights of the distribution. 
This worked relatively well for the ID BHM, and in partic- 
ular for the spectrum-integrated quantities but it may not be 
suited for all possible kind of quenches. We notice that in the 
case of quenches with a mean energy deep in the bulk of the 
spectrum, a generalization of the Lanczos algorithm |76] that 
works in the bulk of a spectrum could be used to get the main 
weights. For what we call small quenches in the following, 
the larger weights are in the low-energy region so Lanczos 
can generically give good results in such situations. 



II. ARGUMENTS ON FINITE SIZE EFFECTS AND THE 
DIFFERENT REGIMES OF A QUANTUM QUENCH 

The goal of this part is to quantify the distance of the 
quench distribution p{E) from the many-body ground-state 
and the low-energy region of the spectrum. A first distance is 
defined from an energetic argument and a second one from the 
overlap with the ground-state of "H/. Both criteria leads to a 
crossover number of bosons Nc{X) that can be computed nu- 
merically and that diverge with small A as a power-law. When 
A^ ^ Nc, the quench probes the low-energy part of the spec- 
trum while when A^ ^ Nc, high-energy physics govern the 
time-evolution. Both definitions do not depend on the integra- 
bility of the model but we may argue that for non-integrable 
models, there is a strong qualitative difference between the 
low-energy part of the spectrum and the bulk of the spectrum. 
These finite-size effects are rather generic while other kind of 
finite-size effects can emerge for a given model: this will be 
for instance the case for the BHM at large U. 
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A. Crossover number of particles from an energetic argument 

The low-energy part of the spectrum - We first have to spec- 
ify what we mean by the low-energy region of the spectrum: 
it corresponds to the typical energies of a few elementary ex- 
citations above the ground-state. These elementary excita- 
tions are quasi-particles, collectives modes, particle-hole ex- 
citations. . . Single or few excitations give a structure (disper- 
sion relations, continuum of low-lying excitations) to the low- 
energy part of the many-body spectrum (see an example in 
Fig.[T]l- We denote by A/ the typical energy scale of a single 
excitation, it is a microscopic energy scale. In Bethe-ansatz 
solvable or free systems, a high energy excitation can be un- 
derstood as a superposition of single-particle excitations while 
this is no longer true for non-integrable systems ll69l pTdh . If 
the number of elementary excitations remains small enough, 
they may hardly interact and have integrable-like features in 
the low-energy part of the spectrum. We thus expect a smooth 
crossover between integrable-like and non-integrable like be- 
haviors with increasing the energy above the ground-state, but 
the typical energy of this crossover is hard to evaluate, except 
that it must be above i?o,/ + A/. 

Criteria - We consider that the energy distribution p{E) 
is centered around the the mean energy E = (V'o.il^/lV'o.i) 
of the distribution (fixed by the initial state) as in general 
AE/{E~Eaj) ^ 1/\/N. Since is not an eigenstate of 
T-Lf, we necessarily have E > Eqj. The criteria we choose to 
distinguish between low-energy (or small) quenches and high- 
energy (or large) quenches is E = E'j (see Fig.|2]i where E^ 
is such that Ej — Eqj = A/ with the ground-state energy 
Eqj. It corresponds to the situation where the mean-energy 
put into the system excites roughly only one elementary exci- 
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FIG. 1: A typical many-body spectrum of a finite-size system: this 
example is taken from the ID Bose-Hubbard model with Uf/J — 
2.5 and L = A'^ = 10. Energies are given as a function of the total 
momentum k. The width of the spectrum is typically proportional to 
or N'^ depending on the statistics. Zooms in the low-energy region 
and in the bulk of the spectrum (grey region) are given. The low- 
energy region features elementary excitations up to a typical energy 
scale Af which is assumed to be microscopic, i.e. non-extensive. 
Here, we take Af — Uf and the relation dispersion of the excitation 
branch is sketch (the line is a guide to the eyes). 




Energy 



FIG. 2: (Color online) Sketch of the energy scales in a quantum 
quench. The initial state builds up an energy distribution p{E) (di- 
agonal ensemble) around a mean energy E fixed by the initial state. 
The quench amplitude A tunes both E and the ground-state energy 



Eoj and the low-energy scale E^. Af = Ej — Eoj is assumed to 
be non-extensive while Eqj and E are extensive. E = Ef defines 
the crossover number of particles Nc- In the thermodynamical limit, 
one expects E ^ EJ for any finite A. 



tation and is thus a finite-size effect. Another way to introduce 
the same criteria is the following: {E — i?o,/)/A/ is the en- 
ergy difference between the initial state and the final ground- 
state in units of the typical elementary excitation energy A/, 
the criteria corresponding to a distance of one A/ |87]. The 
criteria thus amounts to a lower bond of the energies at which 
one enters in the bulk of the spectrum. The order of magni- 
tude of A / is set by the microscopic units of energy of the 
model. For instance, we'll take Uf in the BHM as it controls 
the sound velocity in the superfluid region and the Mott gap 
in the Mott phase. 

This criteria gives a relation between the crossover num- 
ber of particles Nc (on a lattice we work at finite density so 
it also corresponds to a crossover length Lc) and the quench 
amplitude A is such that, if iV ^ Nc{X), the energy is 
mostly distributed among the low-energy excitations while if 
N ^ Nc{X), most of the weights are on high-energy exci- 
tations. We can rewrite the criteria in a more tractable way: 
using the notation e = E/N for the energy per particle, and 
the label for ground-state energies, it reads 



A/ (A) 



e(A) - eo,/(A) 



(5) 



Interestingly, we expect A^c(A) to generically diverge as A~^ 
in the limit of small A. Indeed, we have e = eqj + Xhu 
with hi^i — {ipQ.ilHilipQ.i) / N , the expansion cqj ~ eo.i + 
{deQ/dX)iX + {(PeQ/dX^)iX^ /2, and {deQ/dX)i = hi after 
Feynman-Hellman theorem. With Eq. ([T]), one finally gets 
Nc{X)X^ ^ 2Ai/(d2eo/dA2)j at A ^ 0. 
A few comments can be made on the criteria: 

• When comparing quenches from the same initial state 
but with different "H/, A controls the mean energy per 
particle put into the system. Thus, A is a meaningful 
parameter even in the thermodynamical limit. 

• This definition looks qualitative due to the rather arbi- 
trary choice of A/ and to the fact that, on finite systems, 
the energy distribution can have a rather large width as- 
sociated with energy fluctuations AE. We point out that 
Nc is a crossover number so that N c:^ Nc has no par- 
ticular meaning. Furthermore, from the divergence at 
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small A, one can have 1 <ti N <^ Nc, i.e. a situation 
where energy fluctuations vanish. 

• When A is scanned from to a finite value, both the 
mean energy and the region of the spectrum that plays 
a role in the time-evolution (around e) are continuously 
changed. One can also notice that a quench that starts 
from a ground-state does not necessarily allow to access 
any energy of the H j spectrum, contrary to the situation 
where one prepares the initial state at will. 

• The regimes N ^ Nc and N Nc are expected to 
be physically different for generic (non-integrable) sys- 
tems. Below Af, the density of states is usually much 
smaller than in the bulk of the spectrum: level spac- 
ings are of order of and observables can strongly 
fluctuate with the eigenstate number as it can be seen 
in Figs. |5] and |6] (similar observations can be made in 
the figures of Refs. [23, [12, HH). In this low-energy re- 
gion, RMT arguments are not expected to work (6^ and 
the eigenstates may not be "typical" so we expect the 
ETH to fail. These qualitative observations support the 
difference between the low-energy region and the high- 
energy region of the spectrum made at the beginning of 
this section. 

As the full spectrum width grows as N or N'^ (depend- 
ing on the statistics of the particles) while the number of 
eigenstates grows exponentially with N, the density of 
states in the bulk of the spectrum is exponentially large. 
In this "high-energy" regime (with respect to elemen- 
tary excitations), semi-classical and RMT arguments 
are beheved to work reasonably well for non-integrable 
models ll69ll . which was checked on some strongly cor- 
related svstems i70jl . A s observed numerically on sev- 
eral examples ll27l l32l [33ll . simple observables hardly 
depend on the eigenstate number in this regime, sup- 
porting the ETH. 

• In the thermodynamical limit, we always have N ^ Nc 
and the small quench regime is thus expected to vanish. 
If one wants to check the ETH on a finite size systems, 
one needs sufficiently large A in order to try to reach 
the bulk of the spectrum. However, we will see in this 
paper a counter-example (the BHM) where ETH fails at 
large A (see also Ref.pSh. Even though it looks difficult 
to use quenches to probe very low-energy excitations in 
a very large system, on a finite system, one could tune 
the mean energy from the low to high energy part of 
the spectrum using A. Furthermore, this small quench 
regime is certainly of interest for numerical simulations, 
and also for experiments using a relatively small num- 
ber of atoms (few hundreds or thousands). 

• Lastly, it could be interesting to compare this criteria 
with the domain of validity of bosonization 10, ll^l and 
conformal field theory 10] but this is beyond the 
scope of this paper We note that conformal field theory 
can describe accurately quenches in certain integrable 
models in the thermodynamical limit and for arbitrary 
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FIG. 3: (Color online) Crossover number of bosons Nc obtained 
from the energetic argument in the ID BHM at filling n = 1 and 
starting from Ui = 2. A few points obtained from the criteria based 
on the static fidelity are also given. 



quench ampUtudes lol fToll . Non-integrable models low- 
energy features that are described in terms of a free par- 
ticles (integrable) theory, as bosonization, should dis- 
play non-thermaUzed features as for integrable models. 
In this respect, Ref.'59' gives interesting examples on the 
applicability of these methods to the quench situation. 

We now give examples of Nc{X) for the two models under 
study. In the free boson model, the mean energy after the 
quench can be computed analytically: 
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with eQ i/f = hw^i j 12. The energy fluctuations are given by 
Ae = (e — eQ^i)y'2/N , showing that the distribution gets 
peaked in the thermodynamical limit with the usual scaling. 
A natural choice for A/ is hwf (the only microscopic energy 
scale) and the crossover number of bosons can be expressed 
as a function of the quench amplitude: 



Nc 



hid 



f 



eo,/ 



UJi 



^A + 1 
A2 



This expression diverges as A/X^ in the small quench regime 
and vanishes as 4/ A in the large quench regime. 

For the ID BHM, we take Af = Uf and Nc is given in 
Fig. E] for the particular initial value Ui — 2. It displays the 
expected A^^ divergence at small quenches. We notice that 
the finite size effects on this energy-based criteria are pretty 
small. This can be put on general grounds for ID systems: 
for critical systems the finite-size effects on the ground-state 
energy per particles have a universal correction OTTIl : 



eo(i) = eo(oo)- — + 1 — 



with u the sound velocity and c the central charge. If the sys- 
tem is gapped, the corrections are even smaller as they are 
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exponentially suppressed, by a factor exp(— L/^) with ^ the 
correlation length enters in the formula. In the large quench 
limit of the BHM, one can argue that Nc saturates to a fi- 
nite value. Indeed, in the limit of large A, one finds that 
iVe ^ 2/((n2)o,, - (n2)o,/) + 0(1/A) ~ 2/{n^)o,^, as the 
density fluctuations {n'^)oj are suppressed in the Mott phase. 
Notice that the energy fluctuations, that scale as N~^/'^ in the 
ID BHM, have been computed numerically in Ref. 0. The 
full curve and the two asymptotic behaviors can be simply 
computed from ground-state calculations. 



B. Crossover number of particles based on the static fidelity 

In the thermodynamical limit, the (squared) fidelity be- 
tween the two ground-states F = |(?Ao,i|V'o,/)P is generally 
expected to vanish exponentially with the system size or num- 
ber of particles. Interestingly, l — T counts the contribution of 
the excited states to the time-evolution. A possible definition 
of a crossover number of particles can thus be the value of A 
and N such that T = 1/2, i.e. half of the total weight in the 
ground-state and half in the excited states. In the limit A — > 0, 
one can introduce the fidelity susceptibility Xi.L through the 
expansion ~ 1 — X^Xi,L/2. The scaling of Xi,L is in gen- 
eral non-trivial. If Hi is gapped, the scaling Xi.L ^ L has 
been proposed |78], which gives the divergence Nc ^ A^^. 
In critical systems, super-extensivity, corresponding to a scal- 
ing Xi,L/L ^ i"' with ai > 0, can occur iFtsIItqIi . leading 
to a slower divergence Nc ~ that depends on the 

initial state. Notice that we qualitatively expect that the Nc 
from the fidelity will be smaller than the one based on ener- 
getic argument because, on sufficiently large systems, T can 
be very small while the mean energy is still in the low-energy 
part of the spectrum. 

For the free boson model, the static fidelity as a function of 
A is J- = {VlTX/{l + A/2))^. Setting J" = 1/2, one has 
the crossover number of bosons Nc'. 
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Notice that it also diverges in the small quench regime as 
Nc = 8 In 2/A^ with the same power-law as for the energetic 
arguments. Put in other words, this means that the many-body 
ground-state occupation is robust within a 25% change in w 
for N = 10^, 7% for TV = 10^ and 2% for N = 10^ (see 
next section for the single-particle level occupation). In the 
large amplitude limit, it decreases only logarithmically with 
A, A^c — 2 In 2/ hi A but the prefactor is akeady small. 

The fidelity can also be computed for the ID BHM by Lanc- 
zos calculations. Using the curves J-{X) obtained numerically, 
we determined Nc{X) for numbers of bosons from 6 to 15. 
The result is plotted in Fig. [3] Due to the relatively small sizes 
accessible with Lanczos, we cannot investigate the scaling ex- 
ponent of the small quench divergence. The ground-state fi- 
delity of the ID BHM has been studied in Ref. 80. We observe 
that the static fidelity could be computed on larger chains with 



C. Quench and transition temperature to the Bose-condensed 
regime in the free bosons model 



The free bosons model undergoes a transition to a Bose- 
condensed state below a critical temperature Tc- In the ID 
harmonic trap and on a finite size system, the lowest single 
particle level occupation (no) becomes of the order of N be- 
low Tc ~ HluN/ \n{N) (standard calculations of Tc are per- 
formed in the grand-canonical ensemble and one sees that 
for fixed effective density ljN and iV —> cx), Tc —?> in 
agreement with the fact that their is no Bose-condensation 
in this model in the thermodynamical limit although con- 
densed and non-condensed regimes are clearly seen on finite 
systems). This critical temperature corresponds to a criti- 
cal energy Ec — Eq ^ HujN'^. These standard results can 
be used to answer the question : whether or not a large 
quench from the many-body ground-state can drive the sys- 
tem into the non-condensed regime? We found that the mean- 
energy put into the system scales as E ^ Eqj + hojfNX so 
that A ^ iV is required to reach Ec and the non-condensed 
regime. This surprising behavior (diverging with the num- 
ber of bosons) actually agrees with the exact scaling of the 
single-particle ground-state occupation number which can be 
computed for the quench since we have seen that the distribu- 
tion is the multinomial one : we have (no) ~ Npo ^ N/^/X 
at large A. Similarly, the fluctuations can be computed and 
read (nQ — (no)^) = Npo{l — po) so that the relative fluctu- 
ations scale as 1/ ^/N with a A-dependent prefactor. Conse- 
quently, starting from the many-body ground-state (for which 
(?io) = N), one stays in the condensed regime for finite A and 
one needs A ^ N' with z > 2 to make (no) scale to zero 
in the thermodynamical limit. The physical origin of the fact 
that the quench process makes it difficult to reach the critical 
temperature is that the many-body ground-state has vanishing 
overlaps with the excited states above Tc because they have 
negligible contributions from the single-particle ground-state. 
Starting from a finite temperature state, the quench could help 
cross the critical temperature. 



m. DIAGONAL ENSEMBLE AND THERMALIZATION 

In this section, we compare averages of the expectation val- 
ues of observables obtained from different ensembles: the di- 
agonal, microcanonical and canonical ones. We also show the 
behavior of some local and global observables as a function of 
the energy per particle to discuss the possibility of thermaliza- 
tion according to the ETH. The first numerical evidences that 
the ETH does not work for large quenches on finite systems 
of the ID BHM were recently given in Ref.lssl 
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FIG. 4: Logarithm of the density of states g{e) as a function of the 
energy per particle e in the ID BHM with density n = 1 for three 
different interactions. Mott gaps develop at large U, splitting the 
density of states into many lobes separated by 1 /TV. 



A. Microcanonical temperature and the density of states 

As a preliminary, we discuss the finite size effects and pos- 
sible issues with the microcanonical ensemble in the model 
under study. The standard way to define the microcanonical 
temperature Tm of a closed system is from Boltzmann's for- 
mula 



1 



dsM 
de 



(7) 



where we use the entropy per particle sm — Sm/N and the 
statistical entropy Sm{E) = ki, In n{E). il{E) is the number 
of states within a small energy window 6E around E. Any 
distribution that is peaked enough (6E/E — > in the thermo- 
dynamical limit) will pick up the local density of states g{e) 
through fl{E) ~ g{E)SE. Usually, SE is taken as the en- 
ergy fluctuations with SE ^ E/^/N. Thus, SE is typically 
much larger than microscopic energy scales such as Af. For 
the free boson model, energies per particle are separated by 
huif/N and the degeneracy g{e) of each level can be com- 
puted numerically for small systems. Asymptotic analytical 
results exist in the large energy limit for g{e) llO - ISSll . We can 
thus have access to the microcanonical entropy per particle 
through Sm — ln3(e)/-^- 

In Fig. m we show the logarithm of the density of states of 
the ID Bose-Hubbard model on a finite size system (L = N = 
10) for increasing values of the interaction C/ as a function of 
the energy per particle in units of U. For small interactions, 
the behavior is smooth and one may safely take the derivative 
to get the microcanonical temperature. The system has a den- 
sity of states typical of a bound spectrum Hamiltonian, dis- 
playing first positive and then negative temperature regimes. 
For U = 12 J, in the Mott phase, one observes a gap to the 
ground-state in the low-energy part of the spectrum and also 
some oscillations over a typical scale 1/A^. These oscillations 



are easily understood in the atomic limit (J = 0) where they 
correspond to Mott peaks that have a high degeneracy, giving 
this macroscopic density of states at the center of the lobes. A 
small J broadens the peaks but the lobes are expected to sur- 
vive for large enough [/ in a finite system, as one can see for 
U = 20 J. In this large-?7 limit, eo/U gets close to zero while 
the maximum energy per site is proportional to the number of 
particles (in Fig.lH the situation at high energies is a bit differ- 
ent because we cut the maximum number of bosons onsite). 
The number of Mott lobes being of order iV^, the density of 
lobes per unit of e/U grows as (this remark remains valid 
with a cut in the maximum number of bosons per site). This 
means that the density of states, as a function of the energy 
per site, will be a curve carved into more and more lobes as N 
increases. For large enough systems, Se will be much larger 
than the inter-lobe distance and will pick up the envelope of 
the lobes as a local density of states. On finite systems, Se and 
1 /N could be of the same order of magnitude, which makes 
the definition of the microcanonical temperature rather diffi- 
cult since it is very sensitive to the choice of Se and the shape 
of the peaked distribution. 

In the following, the microcanonical ensemble density- 
matrix pm is defined in the usual way: 



-SE] 



(8) 



with the "free" parameter SE as a "small" energy window en- 
ergy. f2 is simply the number of eigenstates in the energy 
window [E — SE, E + SE]. The sum over the eigenstates of 
H f must be taken over all symmetry sectors. Notice that SE 
can be chosen by hand |27, 32, 33], or in the same way as the 
effective canonical temperature will be: by looking for an ap- 
proximate solution of the equation E = TiIpmH /] (we recall 
that E = (V'o,j:|^/|V-'o,i) is fixed by the initial state). In that 
case, the solution can be multi-valued so it does not necessar- 
ily help. Taking SE as the computed energy fluctuations does 
not help either because on finite systems, the distributions for 
the ID BHM are quite asymmetric and have large moments. 
The choice of SE is in general arbitrary and we have tried to 
choose the one that gives best results for both the correlations 
and the energy. A partial conclusion is that number of parti- 
cles required to have a reliable definition of the microcanon- 
ical ensemble can vary a lot depending on the model and the 
chosen parameters. For the ID BHM, we see that the pecu- 
liar shape of the density of states can be an issue, although it 
intimately linked to the physics of the model. 



B. Canonical ensemble and effective temperature 

Even though we work on a closed system, we introduce a 
canonical density-matrix as it was done in Ref . l30l l32[ l33l and 
implicitly in the (grand)-canonical calculations of Ref. [lit 
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with Z = Tr[e-^^/'=«^«] (9) 



The effective canonical temperature Tb can be defined, as 
in Refs. [3^ [12, [HI as the solution of the equation E = 
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FIG. 5: Local observable 171(e) of the ID BHM as a function of the 
energy per particle for TV = L = 10 and increasing interactions (for 
large quenches, the results were first given in Ref.lSSh 



TtIpbH /] • As the mean energy is a continuous and increasing 
function of Tb, the solution is unique and the optimization 
procedure works well. We take fc^ — 1 in the following so 
that temperatures are given in the same units as the energies. 
Here again, the trace is taken over all symmetry sectors. The 
diagonal ensemble, on the contrary, has non-zero weights only 
in the initial state symmetry sector, that is the even parity sec- 
tor for the free boson model and the fc = sector in the ID 
BHM. As the clouds of points of the distributions sometimes 
look exponential, another temperature can be defined by fit- 
ting the cloud of data with a normalized Boltzmann law and 
using a procedure that minimizes the following cost function 
between two distributions pi and p2'. 



X{Pl,P2) 



Once convergence is reached, we call Tb the effective tem- 
perature obtained from the distribution. 

We lastly recall that provided the density of states scales 
exponentially with the energy and the energy fluctuations are 
negligible in the thermodynamical limit, the microcanonical 
and canonical ensembles will lead to the same thermodynamic 
functions, and same temperatures. 
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FIG. 6: Global observable nk=o{e) as a function of the energy per 
particle (same parameters as in Fig.jSjl. 



where \ipf{e)) is the eigenstate of energy e. gr{e) is a local 
observable since, for a given r, it can be attributed to a sub- 
system. On the contrary, the momentum distribution rik inte- 
grates information from all distances and may be considered 
as a global quantity: 



nfc(e) 



E 

— L + l 



Akr 



9r{e) 



(11) 



In Fig. |5] and |6] one observes that both gi (e) and n^^o (e) 
evolve smoothly in the superfluid regime {U / J = 2.5). One 
also realizes that the largest fluctuations are found in the low- 
energy part of the spectrum, supporting the energetic argu- 
ment for the finite size effects. If one were able to choose e 
in the bulk of the "superfluid" spectrum, one would possibly 
find agreement with ETH. However, for the finite size sys- 
tems at hand, one cannot reach the bulk of the spectrum be- 
fore the Mott lobes emerge with A. As it was shown in Ref.lssI 
and here confirmed, the observables strongly vary within each 
Mott lobe. We now turn the nature of the distributions for dif- 
ferent quenches and compare the results for obtained by the 
different ensembles. Fig. [T] and [8] gather the data for a small 
and large quench from the superfluid region with Ui — 2. 



C. Comparison of observables from different ensembles 

We here focus on the comparison of observables obtained 
from different ensembles in the ID BHM. The evolution of 
one local and one global observable as a function of the eigen- 
states energy per particle is given in Fig. |5] and |6] Each of 
these two observables are used separately in the literature so 
we here give results for both for completeness. The observ- 
ables are the one-particle density-matrix, defined for a trans- 
lationally invariant Hamiltonian as: 

1 ^ 

9r{e) ^ 2^{M^)\bUMi^fie)) , (10) 
1=1 



1. Small quench regime in the ID BHM 

When Uf — 2.5, the distribution is peaked on the final 
ground-state with a large weight po- The tail displays an 
exponential-like behavior that, however, has an effective tem- 
perature Tb different from Tb, determined from the energy. 
This is easily understood from the fact that only the very few 
first weights significantly contribute to the energy, and they 
are not aligned with the tail. As e is very close to e^j in 
this regime and as there are only a very few energies at the 
bottom of the spectrum, the microcanonical ensemble gives a 
bad mean energy and has only a few number of eigenstates. 
In this regime where po is close to one, a minimal micro- 
canonical ensemble would simply be |^o,/) (^o,/|> although 
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FIG. 7: (Color online) Comparison of different ensembles for dif- 
ferent quenches. The effective temperature Tb (given in Fig. [8j is 
fixed by the mean energy. The results are obtained on a system with 
L = N =10. 



it has no statistical meaning. Looking at the correlations g,. in 
Fig.[8]shows that they seem to be thermalized in the sense that 
Pb gives a reasonable account of the correlations. However, 
\ipoj) (V'0,/1 also gives a reasonable account for the correla- 
tions while pm does not satisfactorily reproduce them. The 
system is in a regime dominated by finite size effects, far be- 
low the crossover number of bosons. The points of Ref. [Tl] 
in the "thermalized" region of the phase diagram seem to be- 
long to this regime dominated by finite size effects. We have 
also looked at a slightly larger quench amplitude with Ui ~ 1 
and [// = 4 as in Fig. 3 of Ref. |Tl] (however, we work on a 
slightly smaller system size and the data displayed in Ref. 
were averaged over time, so correlations cannot be quantita- 
tively compared). Since po is smaller, there is a substantial 
difference between the correlations in the final ground-state 
U f and the one from the diagonal ensemble. The canonical 
ensemble still gives the best agreement with p. In a sense, 
the shape of the distributions as given in Ref. 54 does explain 
the observation of Ref. [ill. Yet, the distribution is clearly not 
a true Boltzmann one as the temperature obtained from the 
mean energy and other observables are not identical. In order 
to investigate this deviation, or difficulty to define an effective 
temperature, we have computed the ratio between the two ef- 
fective temperatures Td and Tb in Fig. |9] For L = 6 to 9, 
it remains between 1 and 3.5 and has a tendency to diverge 
at small quenches. Consequently, ETH does not apply here 
due to the presence of strong finite size effects, but one cannot 
claim either that the system is thermalized even though some 
correlations look thermalized in the canonical ensemble. The 
observed distributions are specific to this model and to these 
system lengths and parameters. We also point out that a sim- 
ilar regime has been observed in Ref. 32, corresponding to 
low effective temperatures, but for which the diagonal ensem- 
ble distributions were not plotted. Still, the behavior of large 
systems {N > Nc) in the small quench regime remains an 
open but very interesting question as the low-energy physics 
will control the behavior In this respect, we draw an argu- 



FIG. 8: (Color online) Comparison of averaged correlations Qr cor- 
responding the parameters of Fig.|7] 




FIG. 9: (Color online) Ratio of effective canonical temperatures ob- 
tained from the distribution (To) and from the mean energy (Tb)- 
Results are obtained for a quench starting at Ui = 2. 



ment in favor of non-thermalization: for symmetry reasons, 
the quench only excites states in the ground-state symmetry 
sector while the statistical ensembles average over all symme- 
try sectors. For instance, a system with a branch of excitation 
E{k) can have a A; = gap while the whole spectrum is gap- 
less, hence it could not look thermalized. Starting from a finite 
temperature state or including symmetry breaking terms, like 
disorder, could partially cure this symmetry constraint. 



2. Large quench regime 

Results for two large quenches at a commensurate density 
n — 1, from the superfluid parameters to deep into the Mott 
limit and reversely, are given in Fig. |7] and Fig. [8] For the first 
one, from Ui ~ 2 to U f = 20, the distribution shows very 
strong fluctuations of the weights within each Mott lobes |54]. 
In particular, large weights are present in the low-energy part 
of the first sub-bands. In Ref. 05, it was shown that the larger 
values of gi were correlated to the larger weights (see an- 
other example of such a plot for an incommensurate density 
in Fig. [TTT i. explaining that the ETH does not apply in these 



11 




0.25 0.5 0.75 1 1.25 



^ - 



\\\S\\\\A 








FIG. 10: (Color online) Quench from C/i = 2 to [// = 20 for an 
incommensurate density n = 2/3 for which there is no equilibrium 
quantum critical point. The structure of the density of states, the evo- 
lution of the local correlations gr and the shape of the distributions 
are very similar to the commensurate case. 



finite size systems. This is confirmed by looking at the time- 
averaged correlations that are neither reproduced by p^/ nor 
hy pb- DMRG calculations 1 17, "551] gave evidence that a non- 
thermal correlations gr survive for system sizes of order 100. 

We now elucidate the origin of the observed non- 
thermalized regime, first by looking at the effect of the com- 
mensurability of the density in order to determine whether the 
presence of an equilibrium critical point plays a role for large 
quenches. As shown in Fig.[TO]and[TTl the phenomenology is 
very similar to the commensurate case with a non-thermalized 
regime at large quenches, except that there is no gap above the 
ground-state. Quenches that remain in the superfluid region 
(data not shown) also have the same behavior as for the com- 
mensurate case. These results suggest that the reason for non- 
thermalization is not related to the features of the low-energy 
spectrum, i.e. to the presence of a gap above the ground-state, 
but is related to the proximity of the U — oo limit of the 
model. However, in the small quench regime where the low- 
energy part of the spectrum governs the out-of-equilibrium 
physics, the opening of a gap can certainly play a role. Unfor- 
tunately, due to the finite size effects discussed in this paper, 
this interesting question cannot be addressed with reliability. 
For instance, it has been shown recently [61i] that a quench in 
the quantum Ising model, which is integrable, is sensitive to 
the presence of the critical point. We note that the lobes could 
be qualitatively interpreted as stemming from a ID gapped 
single-particle dispersion relation both in the commensurate 
and incommensurate regimes. However, in the latter case, 
there will not be any transition to an insulating state as a func- 
tion of temperature. 

One can actually argue that the large-?/ structure of the dis- 
tribution is reminiscent of the atomic limit [/ = oo in which 
we show that both the weights and the observables fluctuate 
and are correlated so that ETH is violated in this limit. What 



one can show is that the weights of a quench from /7i = to 
Uf = oo depend on the configuration in each of the degen- 
erated Mott peaks of the Uf — oo limit. This argument does 
not rely on the n = 1 commensurability condition. Indeed, 
the eigenstates of the final Hamiltonian are simply the set of 
configurations ^ with rtj the onsite occupations. The 

energy per particle of the configuration is 
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The initial ground-state is the superfluid state that has equal 
single-particle probabilities pj = 1/L on each site. Using 
formula (|4]i, we get for the diagonal weights: 
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Pn 



(12) 



This makes a connection to the free boson model that we also 
study, having the Uf energy spacing between the degenerate 
levels instead of fkuf and a different energy-configuration re- 
lation. The formula is valid for bare configurations, i.e. when 
they are not symmetrized. Using symmetries, formula (fT2l) 
picks up an additional factor depending on the degeneracy of 
the generalized Bloch state. One can see by taking an exam- 
ple of two configurations with the same energy, or check nu- 
merically, that the weights can be different for configurations 
with the same energy, in the same way as for the free boson 
model. Consequently, in a strongly degenerate Mott peak, the 
diagonal weights are not equal and fluctuate. As soon as a 
non-integrable perturbation (here the hopping J) is turned on 
and lifts the degeneracy, the distribution of the weights will 
still strongly fluctuate within the Mott lobe. This explains the 
findings of Refs.'54','55'and of Fig.]?] Another simple observa- 
tion in this limit is that two degenerate configurations can have 
different expectation values for the observables. An obvious 
one is the onsite particle distribution that counts empty, single, 
double occupations and so on. The off-diagonal correlation 
can be non-zero if the configurations are symmetrized and one 
can check numerically that they actually strongly differ for de- 
generate states. Notice that, in principle, one has to take into 
account the off-diagonal part of the time-averaged density- 
matrix that is non-zero in this highly degenerate limit. When 
one turns on J, this off-diagonal part vanishes and the gr still 
fluctuate strongly for eigenstates close in energy. Lastly, the 
asymmetrical correlation between the weights p„ and the ob- 
servables is also observed in this limit. We show this numer- 
ically on a system with Ui — and Uf — 100 in Fig. [TT| 
(we take Uf/J= 100 and not Jf — because one needs a 
finite, yet very small, J to make p diagonal). The numerics 
for a small J/U in Fig. Q and Fig. |5] strongly supports this 
mechanism as an explanation for the behavior of both the dis- 
tributions and the observables. We remark that the argument 
works as well for the 2D version of the model that was shown 
to have a non-thermalized regime too ifTlll . The fate of this 
explanation in the thermodynamical limit is yet an open ques- 
tion. A scenario could be that this mechanism works above 
a certain critical quench amplitude Xc{N) but how this crit- 
ical value behaves as — > cx) remains a difficult question. 
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FIG. 1 1 : (Color online) Upper panel: comparison of the different 
observables in different ensembles (same parameters as in Fig.llOt. 
Middle panel: the p„ vs gi(n) curve gives proof of the non- 
relaxation towards a thermal state for the same parameters as in the 
upper panel. Lower panel: same plot but for the "integrable" quench 
limit Ui^Q and Uf = 100. 

Consequently, one may understand the finite size effects stem- 
ming from the large-C/ limit as another Nc{\) line in Fig. [3] 
increasing with A and that is specific to this model. Yet, non- 
thermalization in the thermodynamical limit in the BHM can- 
not be excluded as well. Experiments in cold atoms 1 1 ] work 
with a relatively small number of atoms and can easily reach 
this large-[/ limit so that such considerations are physically 
relevant. 

We also give results for a quench from the Mott to the su- 
perfluid limit. There, one could expect from Fig.|5]and |6]that 
ETH could work since the observables behave smoothly with 
e in the final Hamiltonian. However, for the accessible sizes, 
one observes that the Boltzmann law still work better than the 
microcanonical ensemble, with large weights at low energies. 
We conclude that the breakdown of the ETH could here be 
attributed to finite size effects. 



3. Free boson model 

We now briefly discuss the evolution of the distribution 
for the free boson model for a fixed number of bosons 
and increasing A. Very surprisingly, the distribution of the 
single-particle weights versus single-particle energies £2,1 — 
2nhujf + hujf /2 has some remarkable features (we recall that 
only the even levels can be occupied from symmetry reasons). 
In the limit of large energy e ^ 2n, we have 

P2n ^ Po(A) OC ^— 



which has an exponential tail with the effective temperature 
T(A) = hujf/\n In the limit of small quenches, the 

distribution is Boltzmann-like with a temperature r(A) ~ 
— /!W//ln|A/2| going to zero. This exponential-like behavior 
is not generic and a simple counter-example can be found in 
the case of an expanding box |74J. For the many -particle sit- 
uation with = 18 and Ng = 22, we give in Fig.[T2]the evo- 
lution of the distribution for increasing A. For small quench, 
the behavior looks like Boltzmann (we do not expect a pure 
exponential law due to the presence of the degeneracy func- 
tion g{e), see below for a quantitative comparison) and it can 
be understood from the fact that the main contribution comes 
from single-boson excitations that have the same weights as 
the single-particle ones. When A is increased, the energy dis- 
tribution gets peaked around a low-energy level and is strongly 
anisotropic with the maximum at a different place from the 
mean energy. This distribution finally develops a high-energy 
tail for large A. One can compute analytically the third mo- 
ment A/3 = Ti {i5{7i — E)^) which is non-zero and scales as 
N showing that the distribution remains anisotropic and that 
the anisotropy (AIs)^^^ /{E — Eqj) decreases as N^^/^. In 
order to compare the distributions from different ensembles, 
we use the von-Neumann entropy of a density-matrix p which 
is defined as SyN{p) = ^Tr[plnp]. Contrary to observables, 
this quantity is more sensitive to the tail of the distribution. 
SvN /N for the Boltzmann and diagonal ensembles are shown 
in Fig. [13] The density-matrix p'^ is a Boltzmann distribution 
but restricted to the even parity levels only. We see that for 
small quenches, s{p'^) and s{p) are very close. The larger en- 
tropy for s (ps ) is simply due to the fact that half of the Hilbert 
space is not accessible to p for symmetry reasons: s{p'g) and 
s{pb) are actually the same up to a factor 2 in the energy. 
Comparing the data to the microcanonical entropy is not rel- 
evant here because of finite size effects (energy discretization 
and small degeneracy of the first levels) for the values of the 
mean energy accessible here. 
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(e-e^p/co^ 

FIG. 13: (Color online) Von Neumann entropy per particle versus 
energy for the free boson model. Dashed lines are data for A'' = 17, 
full lines for iV = 18 (TV, = 22). 

IV. CONCLUSIONS 

The first conclusion we would like to highlight is that, when 
caiTying out numerical simulations on a finite system, one 
has to care both about the quench amplitude and the size of 
the system to see in which region of the spectrum are the 
main weights of the diagonal ensemble distribution. It has 
been shown that, although the low-energy part of the spec- 
trum is the place where the most interesting physics is ex- 
pected, one experiences large finite size effects when explor- 
ing it. A crossover number of particles, distinguishing be- 
tween the small quench regime and the large quench regime, 
can be tentatively defined from energetic considerations or 
from the static fidelity between the ground-states of the ini- 
tial and final Hamiltonians. One advantage is that they can 
be computed numerically with few finite size effects (for the 
energy based criteria) or with ground-state techniques that 
work on larger systems (for both criteria). The numbers have 
been computed for the two models under study. As the sys- 
tem follows a finite size crossover between the two regimes, 
it can actually happen to be difficult for nowadays numer- 
ics to be close enough to the thermodynamical limit, where 
ETH is expected to work generically, even though some exam- 
ples can be found in the literature |27]. This actually is what 
happens for the one-dimensional Bose-Hubbard model as we 
have seen. Hence, the thermalization-like regime in the small 
quench limit deduced from observables comparison and the 



qualitative Boltzmann-like structure of the distribution cannot 
be considered as truly thermalized because of dominant finite 
size effects. Furthermore, sizes accessible with full diagonal- 
ization cannot reach the bulk of the energy spectrum before 
the structure of the spectrum resembles the infinite-C/ atomic 
limit. The free boson model nicely illustrates the crossover 
from a Boltzmann-like distribution, up to phase space con- 
straints, at small quench to a different distribution. We note 
that due to the large density of states and to negligible energy 
fluctations, we may expect the quench, canonical and micro- 
canonical distributions to eventually be equivalent in the ther- 
modynamical limit. However, we have discussed the fact that 
the smaller the mean-energy (or equivalently the temperature), 
the larger are finite-size effects. We do not believe that the 
observed finite-size and canonical-like distributions at small 
quenches are generic (notice that no claim in that direction 
was made in Ref. l54l) and they may better be simply under- 
stood as (counter-)examples. 

The second important conclusion is that we have shown 
that the non-thermalized regime observed on finite systems 
for large quenches in the ID Bose-Hubbard model is actu- 
ally related to the proximity of the [7 = oo atomic limit, 
something that may qualitatively be equivalent to the prox- 
imity of an integrable point. Indeed, this regime does not 
depend on the low-energy features at a commensurate den- 
sity, i.e. to the presence of the superfluid-Mott transition, and 
besides, the structure of the diagonal ensemble stems from 
the U — Q ^ U f = CO quench limit of the Bose-Hubbard 
model. In non-integrable models, the challenging issues on 
what are the features of the small quench regime for very large 
sizes and how the non-thermalized regime neighboring an in- 
tegrable point survive in the thermodynamical regime seem to 
be hardly accessible to current numerical algorithms. 
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